arXiv:1502.07505vl [stat.ME] 26 Feb 2015 


A mixed effect model for bivariate meta-analysis of diagnostic 
test accuracy studies using a copula representation of the 
random effects distribution 

Aristidis K. Nikoloulopoulos* 


Abstract 

Diagnostic test accuracy studies typically report the number of true positives, false positives, 
true negatives and false negatives. There usually exists a negative association between the num¬ 
ber of true positives and true negatives, because studies that adopt less stringent criterion for 
declaring a test positive invoke higher sensitivities and lower specificities. A generalized lin¬ 
ear mixed model (GLMM) is currently recommended to synthesize diagnostic test accuracy 
studies. We propose a copula mixed model for bivariate meta-analysis of diagnostic test ac¬ 
curacy studies. Our general model includes the GLMM as a special case and can also operate 
on the original scale of sensitivity and specificity. Summary receiver operating characteristic 
curves are deduced for the proposed model through quantile regression techniques and differ¬ 
ent characterizations of the bivariate random effects distribution. Our general methodology is 
demonstrated with an extensive simulation study and illustrated by re-analysing the data of two 
published meta-analyses. Our study suggests that there can be an improvement on GLMM in 
fit to data and makes the argument for moving to copula random effects models. Our modelling 
framework is implemented in the package CopulaREMADA within the open source statistical 
environment R. 

Keywords: copula models; diagnostic tests; multivariate meta-analysis; random effects mod¬ 
els; SROC, sensitivity/specificity. 


1 Introduction 

Synthesis of diagnostic test accuracy studies is the most common medical application of multivariate 
meta-analysis [21, 30]. Meta-analysis is broadly defined as the quantitative review of the results of 
related but independent studies [41]. The purpose of a meta-analysis of diagnostic test accuracy 
studies is to combine information over different studies, and provide an integrated analysis that will 
have more statistical power to detect an accurate diagnostic test than an analysis based on a single 
study. Accurate diagnosis plays an important role in the disease control and prevention [29]. 

Diagnostic test accuracy studies observe the result of a gold standard procedure which defines 
the presence or absence of a decease and the result of a diagnostic test. They typically report the 
number of true positives (diseased people correctly diagnosed), false positives (non-diseased people 
incorrectly diagnosed as diseased), true negatives and false negatives. As the sensitivity (proportion 
of those with the disease) and specificity (proportion of those without the disease) are estimated 
from different samples in each study (diseased and non-diseased patients), they can be assumed 
to be independent so that the within-study correlations are set to zero [30]. However, there may 
be a negative between-studies association which should be accounted for. A negative association 
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between these quantities across studies is likely because studies that adopt less stringent criterion 
for declaring a test positive invoke higher sensitivities and lower specificities [21]. 

In situations where studies compare a diagnostic test with its gold standard, heterogeneity arises 
between studies due to the differences in disease prevalence, study design as well as laboratory and 
other characteristics [7]. Because of this heterogeneity, a generalized linear mixed model (GLMM) 
has been recommended in the biostatistics literature [4, 1, 14, 29] to synthesize information. Note in 
passing that it is equivalent with the hierarchical summary receiver operating characteristic model in 
Rutter and Gatsonis [46] for the case without covariates [15, 5]. The GLMM assumes independent 
binomial distributions for the true positives and true negatives, conditional on the latent pair of 
transformed (via a link function) sensitivity and specificity in each study. The random effects (latent 
pair of transformed sensitivity and specificity) are jointly analysed with a bivariate normal (BVN) 
distribution. 

Chu et al. [7] propose an alternative mixed model which operates on the original scale of 
sensitivity and specificity. The random effects follow the bivariate Saimanov’s [47] family of dis¬ 
tributions with beta margins [28]. However, this random effects distribution has a limited range 
of dependence and is inappropriate for general modelling unless the responses are weakly depen¬ 
dent. Hence, this model is too restrictive in the context of diagnostic accuracy studies where strong 
(negative) dependence is likely. 

We propose a copula mixed model as an extension of the GLMM and mixed model in Chu et 
al. [7] by rather using a copula representation of the random effects distribution with normal and 
beta margins, respectively. Copulas are a useful way to model multivariate data as they account for 
the dependence structure and provide a flexible representation of the multivariate distribution. The 
theory and application of copulas have become important in finance, insurance and other areas, in 
order to deal with dependence in the joint tails. Here, we indicate that this can also be important in 
meta-analysis of diagnostic test accuracy studies. Diagnostic test accuracy studies is a prime area 
of application for copula models, as the traditional assumption of multivariate normality is invalid 
in this context. 

A copula approach for meta-analysis of diagnostic accuracy studies was recently proposed by 
Kuss et al. [27] who explored the use of a copula model for observed discrete variables (number 
of true positives and true negatives) which have beta-binomial margins. This model is actually 
an approximation of a copula mixed model with beta margins for the latent pair of sensitivity and 
specificity. Although, this approximation can only be used under the unrealistic case that the number 
of observations in the respective study group of healthy and diseased probands is the same for 
each study. In real data applications, the number of true positives and negatives do not have a 
common support over different studies, hence, one cannot conclude that there is a copula. The 
natural replicability is in the random effects probability for sensitivity and specificity. 

The remainder of the paper proceeds as follows. Section 2 summarizes the standard GLMM 
for synthesis of diagnostic test accuracy studies. Section 3 has a brief overview of relevant copula 
theory and then introduces the copula mixed model for diagnostic test accuracy studies and discusses 
its relationship with existing mixed models. Section 4 discusses suitable parametric families of 
copulas for the copula mixed model, deduces summary receiver operating characteristic curves for 
the proposed model through quantile regression techniques and different characterizations of the 
bivariate random effects distribution, and demonstrates that they can show the effect of different 
model assumptions. Section 5 contains small-sample efficiency calculations to investigate the effect 
of misspecifying the random effects distribution on parameter estimators and standard errors and 
compare the proposed methodology to existing methods. Section 6 summarizes the assessment 
of the proposed models using the Vuong’s statistic [53], which is based on sample difference in 
Kullback-Leibler divergence between two models and can be used to differentiate two parametric 
models which could be non-nested. Section 7 presents applications of our methodology to four data 
frames with diagnostic accuracy data from binary test outcomes. We conclude with some discussion 
in Section 8, followed by a section with the software details and a technical Appendix. 

2 The standard GLMM 

We first introduce the notation used in this paper. The focus is on two-level (within-study and 
between-studies) cluster data. The data are are {yij,nij), i = 1,..., N, j = 1,2, where j is an index 
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for the within study measurements and i is an index for the individual studies. The data, for study 
i, can be summarized in a 2 x 2 table with the number of true positives (ya), true negatives {yi 2 ), 
false negatives (nn — yn), and false positives (nj 2 — 2 / 12 ); see Table 1. 

Table 1: Data from an individual study in a 2 x 2 table. 


Test 

Disease (by gold standard) 


Yes 

No 

Positive 

yn 

ni2 — ya 

Negative 

'iT'ii yn 

ya 

Total 

nn 

na 


The standard two-level model of meta-analysing diagnostic test accuracy studies [4, 15, 1, 14, 
29] lies in the framework of mixed models [8]. The within-study model assumes that the number 
of true positives Yu and true negatives are conditionally independent and binomially distributed 
given X = X, where X = {Xi,X 2 ) denotes the bivariate latent (random) pair of transformed 
sensitivity and specificity. That is 


YillXy = XI 


— X2 


r\j 


r\j 


Binomial nil, I 
Binomial (ni 2 ,1 



( 1 ) 


where /(•) is a link function such as the commonly used logit. The between studies model assumes 
that X is BVN distributed with mean vector /r = (((tti), l{'it 2 ))^ and variance covariance matrix 

^T^VThatis 

\pa1a2 a2 ) 

X~BVN(/r,S). (2) 

The models in (1) and (2) together specify a GLMM with joint likelihood 

N 2 

L{TTi,TT 2 ,( 7 l,a 2 , p) = u m yij]nij,l ^(Xj)^(I) i2{xi,X 2; fx, Yi)dxidx2, 

1=1-^ j=i 


where 

c/(y;n,7r) = - 7r)"“^, y = 0,l,...,n, 0 < vr < 1, 

is the binomial probability mass function (pmf) and 0i2(-; X) is the BVN density with mean 
vector ft and variance covariance matrix S. The parameters tti and tt 2 are those of actual interest 
denoting the meta-analytic parameters for the sensitivity and specificity, respectively, while the 
univariate parameters af and a 2 are of secondary interest denoting the variability between studies. 


3 The copula mixed model for diagnostic test accuracy studies 

In this section, we introduce the copula mixed model for diagnostic test accuracy studies and discuss 
its relationship with existing mixed models. Before that, the first subsection has some background 
on copula models. In Subsection 3.2 and Subsection 3.3 a copula representation of the random 
effects distribution with normal and beta margins respectively is presented. We complete this section 
with details on maximum likelihood estimation. 


3.1 Overview and relevant background for copulas 

A copula is a multivariate cumulative distribution function (cdf) with uniform (7(0,1) margins [22, 
33, 25]. If Fi 2 is a bivariate cdf with univariate margins Fi, F 2 , then Sklar’s [51] theorem implies 
that there is a copula C such that 

Fi2{xi,X2) = c(^Fi{xi),F2{x2)y (3) 
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The copula is unique if Fi, F 2 are continuous, but not if some of the Fj have discrete components. 
If Fi 2 is continuous and (Xi, X 2 ) ~ F 12 , then the unique copula is the distribution of (?7i, U 2 ) = 
(Fi(Xi),F 2 (X 2 )) leading to 

C{ui,U 2 ) = Fi 2 (^F~^ (ui), F~^ {U 2 )^ , 0 <Uj < l,j = 1,2, 

where F~^ are inverse cdfs. In particular, if <I>i 2 (-; p) is the BVN cdf with correlation p and standard 
normal margins, and <1> is the univariate standard normal cdf, then the BVN copula is 

C{ui,U 2 ) = ‘^I 2 (^^~^{ui),^~^{u 2 ); pY 

The power of copulas for dependence modelling is due to the dependence structure being considered 
separate from the univariate margins; see e.g., [22, Section 1.6]. If C{-; 9) is a parametric family of 
copulas and Fj{-] pj) is a parametric model for the jth univariate mai'gin, then 

C (^Fi (xi; r/i), F 2 (x2; r/2); 0 ) 

is a bivariate parametric model with univariate margins Fi,F 2 . For copula models, the variables 
can be continuous or discrete [38]. 


3.2 The copula mixed model for the latent pair of transformed sensitivity and specificity 

Here we generalize the GLMM by proposing a model that links the two random effects using a 
copula function rather than the BVN distribution. 

The within-study model is the same as in the standard GLMM; see (1). The stochastic represen¬ 
tation of the between studies model takes the form 

(ch(Vi;^(7ri),af),cl>(X2;Z(7r2),c72)) ~C(-;0), (4) 

where C(-; 9) is a parametric family of copulas with dependence parameter 9 and <!>(•; p, a"^) is the 
cdf of the N(/r, distribution. The joint density fi 2 ixi,X 2 ) of the transformed latent proportions 
can be derived as a double partial derivative of the cdf in (3) 

9(7 f 4> (xi; Z (tti), o-f), 4> ( 3:2 ; /(7r2), 

fi2{xi,X2-,Tri,TT2,ai,a2,9) = --—-- (5) 

OX 1 OX 2 

= c(^4>(xi; /(vTi), erf), 4 >(x 2; 1{'K2),cfI) ; ^(vri), erf) (/>(3:2; 1{tt2),cfI) , 

where c{ui,U2;9) = d‘^C{ui,U2\9)/duidu2 and 0(-;/r,cr^) is the copula and N(^,cr^) density, 
respectively. The models in (1) and (4) together specify a copula mixed model with joint likelihood 


L(7ri,7r2,ai,cr2,6») = 


^ poo poo ^ 
i=i "'-oo j=l 


( 6 ) 


4>(3:2;?(vr2),o-f);6» )n (p(xj-, l{7rj),aj)dxidx2 

3=1 

n L L yij-,nij,l ^( 4 > ^{uj-,l{7rj),a'j))]c{ui,u2;9)duidu2. 


i=l 


0 30 


3=1 


It is important to note that the copula parameter 0 is a parameter of the random effects model 
and it is separated from the univariate parameters. The univariate parameters vri and 7:2 are those 
of actual interest denoting the meta-analytic parameters for the sensitivity and specificity, while the 
univariate parameters erf and erf are of secondary interest expressing the variability between studies. 
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3.2.1 Relationship with the GLMM 


In this subsection, we show what happens when the bivariate copula is the BVN copula. The 
resulting model is the same as the GLMM. 

The BVN copula density is 


c{ui,U2;p) = 


\/l - P‘ 


: exp 


+ ^2 ~ ‘^p^l^2\ 
, 2^1^ ) 


exp 


zf + z^ 


where Zj = j = 1, 2. Thenfor Uj = ^(vrj), aj) we have zj = {xj — l{'Kj)) /aj, j = 

1, 2. Hence, the joint density in (5) becomes 


/l2(a^l, X 2 ; TTi, 7r2, (7i, (72, p) — 


exp 


(xi - ^(7ri))" 

27r(7i(72y/l - ^ I 2cr^ 

(x2-^(vr2))^ (xi -/(vTi)) (X2 - Z(vr2)) 




+ 


2cj| 


-2p 


(71(72 


which apparently is the BVN density (^i 2 (ici, X 2 ; pL, S). 


3.3 The copula mixed model for the latent pair of sensitivity and specificity 

The within-study model also assumes that the number of true positives Yu and true negatives Yi 2 are 
conditionally independent and binomially distributed given X = x, where X = {Xi,X 2 ) denotes 
the bivariate latent random pair of sensitivity and specificity. That is 

Yii\Xi = xi ~ Binomial (nil, aJi); 

yi2\X2 = X2 ~ Binomial (ni2,3:2). ( 7 ) 

So one does not have to transform the latent sensitivity and specificity and can work on the original 
scale. The Beta(Q;, /3) distribution can be used for the marginal modeling of the latent proportions 
and its density is 

x^ — x^^ ^ 

f{x;a,l 3 ) = - ——— -, 0 <a;<l, a,P> 0 . 

ijyCX^ Pj 

In the sequel we will use the Beta(7r,7) parametrization, where tt = (mean parameter) and 

7 = a^p+i (dispersion parameter). 

The stochastic representation of the between studies model is 

(F(Xi;7ri,7i),F(X2;7r2,72)) ~C7(-;0), (8) 

where C'(-; 6 ) is a pai'ametric family of copulas with dependence parameter 9 and F{-\ tt, 7 ) is the 
cdf of the the Beta( 7 r, 7 ) distribution. The models in (7) and ( 8 ) together specify a copula mixed 
model with joint likelihood 

(-1 H 2 

L(7ri,7r2,7i,72,6l) = n / n g{yij]nij,Xj)c\F{xr,-Ki,-ii),F{x2]T^2,l2)\0j 

i=i “'o -^0 j=l 

2 

^ n fixj]'^j,'lj)dxidx2 (9) 

i=i 

^ rl p1 ^ 

= IT / / W 9{yij'^'^ij^y'~^iuj'^'^j^'yj))c{ui,U2-,6)duidu2. 

i=Po Jo 

As before, the copula parameter 0 is a parameter of the random effects model and it is separated 
from the univariate pai'ameters, the univariate parameters tti and 712 are the meta-analytic parameters 
for the sensitivity and specificity, but, now 71 and 72 express the variability between studies. 
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3.3.1 Relationship with existing models 

Chu et al. [7], instead of using a copula for the random effects distribution or a copula density for 
X in (9), use the Sarmanov’s [47] family of bivariate densities 

h2{xi,X2) = fl{xi)f2{x2){l + 6'V’i(xi)V'2(x2)), 

where fj{-) is the marginal density of Xj, is abounded non-constant function such as fj{x) 
il’j{x)dx = 0, and 1 -|- 0 il;i{xi)Tp2{x2) > 0 for all xi,X 2 - For the Sarmanov’s densities if one uses 
V’j = 1 — 2Fj{xj), j = 1,2, then the Farlie-Gumbel-Morgenstern copula (density) is obtained. 
However in [7], “kernels” of the type i/jjixj) = Xj — E{Xj) are considered as in [28]. The advan¬ 
tage of this choice is that the corresponding likelihood function has a closed form, since the product 
of integrals can be evaluated analytically. The joint likelihood takes the form 


L(7ri,7r2,7i,72,6») 


^ /.l .1 2 2 

n/ n giVij ; nij , Xj)f{xj + on {xj — •Kj)jdxidx2 

i=l -^0 Jo 


N 2 

YlUHvij-,' 

1=1 j=l 


l + ofl 


Vij njjTTj 

- 1 


i=i 7 


where 


fn\ ^( 2 / + ^/7 - 7r,n - y-h (1 - 7 r)(l - 7 ) 77 ) 

/i(y;n,7r,7) = ^- - -r-y = 0,1,..., n, 0 < tt, 7 < 1, 

V2// 7r/7 — TT, (1 — 7 r)(l — 7 )/ 7 j 

is the pmf of a Beta-Binomial(n, vr, 7 ) distribution with mean nvr and variance n 7 r(l — vr) (l -|- (n — 
1 ) 7 ). The disadvantage of this mixed model is that the Sarmanov’s density with beta margins in [28] 
has a limited range of dependence and is inappropriate for general modeling unless the responses 
are weakly dependent. 

Kuss et al. [27] proposed a copula model with beta-binomial margins in this context. This 
model is actually an approximation of the copula mixed model with beta margins for the latent pair 
of sensitivity and specificity in (7) and ( 8 ). They attempt to approximate the likelihood in (9) with 
the likelihood of a copula model for observed discrete variables which have beta-binomial margins. 

The approximation that they suggest is 

N 2 

L( 7 ri, 7 r 2 , 7 i, 72 , 6 ») c(^iT(yii; na, tti, 71 ), iT(yi 2 ; ni 2 , vr 2 , 72 ); h{yij-,nij,TTj,jj), 

i=l j=l 

where is the cdf of the the Beta-Binomial(n, vr, 7 ) distribution. In their approxima¬ 

tion the authors also treat the observed variables which have beta-binomial distributions as being 
continuous, and model them under the theory for copula models with continuous margins. Kuss 
et al. [27], referring to Genest and Neslehova [9], claim that there are problems on applying cop¬ 
ula to discrete data especially in extreme cases with very small numbers of support points for the 
discrete marginal distributions. Genest and Neslehova [9] only warn against estimation for discrete- 
margined copula models using rank-based methods, instead recommending maximum likelihood 
estimation. Essentially, Genest et al. [10] apply copula models to multivariate binary data (the ex¬ 
treme case of discreteness) and call on composite likelihood techniques for estimation. Multivariate 
copulas for discrete response data have been in use for a considerable length of time, e.g., in Joe 
[22], and earlier for some simple copula models. Several examples of copula models for multivari¬ 
ate discrete data can be found in the literature; see e.g., [40] for an application in biostatistics and 
[34] for a survey of copula models and methods for multivariate discrete response data. 
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However, the main problem in [27] is that the approximation (even if treating the observed 
variables which have beta-binomial distributions as being discrete) can only be used under the un¬ 
realistic case that the number of observations in the respective study group of healthy and diseased 
probands riij is the same for each study i. In real data applications, the discrete Yij do not have a 
common support over different studies or i, hence, one cannot conclude that there is a copula for 
{Yii, Yi 2 ) that applies when the Uij vary with different studies i. The natural replicability is in the 
random effects probability for sensitivity and specificity. 


3.4 Maximum likelihood estimation and computational details 

Estimation of the model parameters (vri, 7 r 2 , ai,a2,0) and (tti, 7 r 2 , 71 , 72 , can be approached by 
the standard maximum likelihood (ML) method, by maximizing the logarithm of the joint likelihood 
in ( 6 ) and (9), respectively. The estimated pai'ameters can be obtained by using a quasi-Newton [32] 
method applied to the logarithm of the joint likelihood. This numerical method requires only the 
objective function, i.e., the logarithm of the joint likelihood, while the gradients are computed nu¬ 
merically and the Hessian matrix of the second order derivatives is updated in each iteration. The 
standard errors (SE) of the ML estimates can be also obtained via the gradients and the Hessian 
computed numerically during the maximization process. Assuming that the usual regularity condi¬ 
tions [49] for asymptotic maximum likelihood theory hold for the bivariate model as well as for its 
margins we have that ML estimates are asymptotically normal. Therefore one can build Wald tests 
to statistically judge any effect. 

Lor mixed models of the form with joint likelihood as in ( 6 ) and (9), numerical evaluation of the 
joint pmf is easily done with the following steps: 

1. Calculate Gauss-Legendre quadrature points {uq : q = 1,... ,nq} and weights {wq : q = 

1.. .. , ng} in terms of standard uniform; see e.g., [52]. 

2. Convert from independent uniform random variables {uq^ : qi = 1,..., n^} and {uq^ : q 2 = 

1 .. .., Uq} to dependent uniform random variables {uq^ : gi = 1 ,..., Uq} and {C~^ (^^92 1 ^?!! • 

= 92 = 1, • • •, riq} that have distribution C'(-; 9). The inverse of the conditional distribu¬ 
tion C{v\u-, 9) = dC{u, v] 9)/du corresponding to the copula C'(-; 9) is used to achieve this. 

3. Numerically evaluate the joint pmf, e.g., 

1 2 

'[\g{yj]nj,F~^{uj\T:j,-ij))c{ui,U2\9)duidu2 
i=i 



in a double sum: 

Uq Uq 

( 2/1 \n,F~^ {uq ^; , 7 ^-))p (^^ 2 ; n, {C~^{uq^\uq ^; 9 ); vr^, 7 ^) ^ 

gi=l 92=1 


With Gauss-Legendre quadrature, the same nodes and weights are used for different functions; 
this helps in yielding smooth numerical derivatives for numerical optimization via quasi-Newton 
[32]. Our comparisons show that = 15 is adequate with good precision to at least at four 
decimal places; hence it also provides the advantage of fast implementation. 

To sum up, our mixed effect model for meta-analysis of diagnostic test accuracy studies using 
a copula representation of the random effects distribution with a double integral is straightforward 
computationally. Note in passing that the linear mixed model in [43] can also provide handy com¬ 
putations, but it has limitations due to the use of continuity correction and normal approximation 
[4, 29]. 
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4 Choices of parametric families of copulas 


In our candidate set, families that have different strengths of tail behaviour (see e.g., [17]) are in¬ 
cluded. In the descriptions below, a bivariate copula C is reflection symmetric if its density satisfies 
c{ui,U 2 ) = c(l — til, 1 — U 2 ) for all 0 < ui,U 2 < 1. Otherwise, it is reflection asymmetric often 
with more probability in the joint upper tail or joint lower tail. Upper taii dependence means that 
c(l — n, 1 — tt) = 0 {u~^) as u —)• 0 and iower taii dependence means that c{u,u) = 0{u~^) 
as u ^ 0. If (f/i, 172 ) ~ C for a bivariate copula C, then {1 — Ui,l — U 2 ) ~ Cigo", where 
0*180° (ni, 1 x 2 ) = -f ri 2 — 1 -I- 0(1 — 1 x 1 ,1 — 1 x 2 ) is the survival (or rotated by 180 degrees) copula 
of C*; this “reflection” of each uniform [7(0,1) random variable about 1/2 changes the direction of 
tail asymmetry. 

• Reflection symmetric copulas with tail independence satisfying 0(xx, u) = O(xx^) and 0(1 — 
XX, 1 — xx) = O(xx^) as xt —>• 0, such as the Frank copula with inverse conditional cdf 


O ^(x;|xx;0) 



l-(l-e-®) 
(x;“^ — l)e“®“ + 1 


9 G (—00,00) \ { 0 }. 


• Reflection symmetric copulas with intermediate tail dependence [18] such as the BVN copula, 
which satisfies 0(xx, xx, 9) = (7(xx^/^^+^) (— log xx)“^0i+®)) as u —)■ 0 with inverse conditional 
cdf 


O ^(x;|xx; 0) = ^(u)^, 0 g[—1,1]. 


Reflection asymmetric copulas with lower tail dependence only such as the Clayton copula 
with inverse conditional cdf 

C-\v\u;9) = |(x;-0(i+®) - 1 )^-® + 9 G ( 0 ,oo). 

Reflection asymmetric copulas with upper tail dependence only such as the rotated by 180 
degrees Clayton copula with inverse conditional cdf 


C-i(x;|xx; 0 ) = 1 - [{(1 - t;)-'^/(i+®) - 1}(1 - xx)"® + 1 


-i/e 


9 G ( 0 ,00). 


The Frank and BVN copulas interpolate from the Frechet lower (perfect negative dependence) to 
the Frechet upper (perfect positive dependence) bound, and, thus they are sufficient from bivariate 
studies on diagnostic accuracy where negative dependence between the number of true positives 
and true negatives is expected. The Clayton copula belongs in the Archimedean class of copulas. 
Archimedean copulas, see e.g. [22], have the form, 

C(xxi, XX2 ; 9) = (j) ; 9) + (j)~^{u2 ] 9) ; 9) , ( 10 ) 

where the generator (f){u ; 9) is the Laplace transform (LT) of a univariate family of distributions of 
positive random variables indexed by the parameter 9, such that (/(•) and its inverse have closed 
forms. The Clayton copula interpolates from the independence {9 —)• 0) to the Frechet upper 
(comonotonic copula) bound {9 —>• 00 ). For extension of the Laplace transform for 9 G [—1,0), the 
Clayton family extends to countermonotonicity {9 —1). However this extension is not generally 

useful for applications because the support of (10) is not all of (0,1)^ [22, page 109]. Negative 
dependence in Clayton copulas can be introduced by applying decreasing transformations to the 
“oppositely” ordered variables. If {Ui, U 2 ) ~ C where C* is a copula with positive dependence, one 
could always get some negative dependence, by supposing 6 * 90 ° is the copula of ((7i, 1 — ( 72 ) (rota¬ 
tion by 90 degrees) or 6 * 270 ° the copula of (1 — (7i, U 2 ) (rotation by 270 degrees). So it is worthwhile 
to rotate the Clayton copula by 90 and 270 degrees to model negative dependence. These rotated 
copulas interpolate from the Frechet lower (perfect negative dependence) {9 —)• 00 ) to the indepen¬ 
dence {9 —?■ 0). Negative upper-iower taii dependence means that c(l — u,u) = 0{u~^) as xx —^ 0 
and negative iower-upper taii dependence means that c(u, 1 — xx) = 0(xx"^) as xx —)• 0 [24]. So in 
order to model negative (tail) dependence the choices are: 








• Reflection asymmetric copula family with negative upper-lower tail dependence, such as the 
rotated by 90 degrees Clayton copula with inverse conditional cdf 



9 e (0,oo). 


• Reflection asymmetric copula family with negative lower-upper tail dependence, such as the 
as the rotated by 270 degrees Clayton copula with inverse conditional cdf 


C 9) = 1 — {(1 — — l]u^ -|- 1 ^ 0 £ (0, oo). 


For this paper, the above copula families are sufficient for the applications in Section 7, since 
tail dependence is a property to consider when choosing amongst different families of copulas and 
the concept of upper/lower tail dependence is one way to differentiate families. Nikoloulopoulos 
and Karlis [39] have shown that it is hard to choose a copula with similar properties from real data, 
since copulas with similar (tail) dependence properties provide similar fit. Kuss et al. [27] used, in 
addition to these copulas, the Placket copula. Plackett copula is a reflection symmetric copula [22, 
pages 221-22] with tail independence [33, page 215] (not reflection asymmetric copula as stated in 
[27]) and is not used here since we have included another choice of copulas with similar properties 
i.e., the Frank copula. 

4.1 Summary receiver operating characteristic curves 

Rutter and Gatsonis [46] proposed a hierarchical summary receiver operating characteristic (SROC) 
curve which for some cases is the same with the corresponding GLMM SROC curve [5]. For the 
GLMM model, the model parameters control the shape of the SROC curve. The GLMM SROC 
curve can be obtained through a characterization of the estimated bivariate normal distribution by a 
line [5, 6, 7]. Based on the bivariate normality of the random effects, the expected sensitivity for a 
chosen specificity in the transformed scale is given in a closed form: 


E[Xi\X2 = X 2 ] = [((vTi) - pl{-K2)cri/a2] + pl{x2)cri/a2. 


( 11 ) 


In general, however, E[Xi\X 2 = X 2 ] is not in closed form and thus does not have simple expressions 
in terms of distribution functions and copulas. 

An alternative to the mean for specifying “typical” values of Xi for each value of X 2 is the 
median, which leads to the notion of median regression of Xi on X 2 . For X 2 in range of X 2 , let 
xi := xi{x 2 ) denote a solution to the equation Pr(2fi < xi\X 2 = X 2 ) = 1/2. Then the scatter plot 
of xi{x 2 ) and X 2 is the median regression curve of Xi on A 2 . 

For copula models, median regression curves [33, pages 217-218] can be easily calculated, 
since 


Pr(Xi < X 1 IX 2 = X 2 ) = Pr(C/i < Fi(xi)|C /2 = E 2 {x 2 )) = C{ui\u 2 ) 


Ui = Ei{xi) 

U 2 = E 2 {x 2 ) 


but their shape also depends on the choice of bivariate copulas. Furthermore, as emphasized in [45], 
since there is no unique definition of a SROC curve, it is preferable and will make more sense to 
deduce confidence regions as well. To this end in addition of using just median regression curves 
we will also exploit the use of quantile regression curves with a focus on high {q = 0.99) and 
low quantiles {q = 0.01) which are strongly associated with the upper and lower tail dependence 
imposed from each parametric family of copulas. These can be also seen as confidence regions 
of the median regression SROC curve. Note that Kendall’s tau only accounts for the dependence 
dominated by the middle of the data, and it is expected to be similar amongst different families 
of copulas. However, the tail dependence varies, as explained in Section 4, and is a property to 
differentiate amongst different families of copulas. 

To find the quantile regression curves: 
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1. Set C{ui\u 2 ] 0) = q. 

2. Solve for the quantile regression curve ui := ui{u 2 , Q] 0) = C~^{q\u2]0). 

3. For j = 1,2 replace Uj by Fj{xj]7rj,^j) for beta margins or ^j(^Xj]l{Trj),aj) for normal 
margins. 

4. Plot xi := xi(x 2 ,q) versus X 2 - 

Of course, the quantile regression curve X 2 '■= X 2 {xi,q) of X 2 on Xi is defined similarly. 
However, there is no priori reason to regress xi on X 2 instead of the other way around [1]. In fact, 
if one wants to reserve the nature of a bivariate response instead of a univariate response along 
with a covariate, then a contour graph can be easily plotted. The contour plot can be seen as the 
predictive region (analogously to [43]) of the estimated pair of sensitivity and specificity. However, 
the resulted shape of the prediction region is not depended on the assumption of bivariate normality 
for the random effects. 



Figure 1: Contour plots and quantile regression curves from the copula representation of the random effects distribution 
with normal margins and BVN, Frank, and Clayton by 90 and 270 copulas with the same model parameters {tti = 
0.7,712 = 0.9, (Ti = 2,(T2 = l,r = —0.5}. Red and green lines represent the quantile regression curves xi ~ 
xi(x2,q) and X2 := X 2 (xi,q), respectively; for q = 0.7> solid lines and for q G {0.01,0.99} dotted lines. 

To depict the different shapes of the SROC curves, in Figures 1 and 2 we plot them from the 
copula representation of the random effects distribution with normal and beta margins, respectively, 
and BVN, Frank Clayton by 90 and 270 copulas with the same model parameters {tti = 0.7, 7 r 2 = 
0.9,fJi = 2,(72 = l,r = —0.5} and {tti = 0.7, 7 r 2 = 0.9 ,71 = 0 . 2,72 = 0.1,r = —0.5}, 
respectively. We convert from r to the BVN, Frank and rotated Clayton copula parameter 9 via the 
relations 

2 

r = — arcsin(0), (12) 

TT 
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r = 


(13) 


1 - 40-1 - 40-2/e , 0<0 

1 _ 40-1 + 40-2 j-e , 0>O ’ 

( 0/(0 + 2) , by 0 or 180 degrees 

\ —0/(0+ 2) , by 90 or 270 degrees 

in [19], [11], and [12] respectively. 


BVN 


Frank 




Clayton by 90 


Clayton by 270 



(14) 


Figure 2: Contour plots and quantile regression curves from the copula representation of the random effects distribution 
with beta margins and BVN, Frank, and Clayton by 90 and 270 copulas with the same model parameters {tti = 0.7, ii 2 = 
0.9, 71 = 0.2, 72 = 0.1, r = —0.5}. Red and green lines represent the quantile regression curves xi := xi {x 2 , q) and 
X 2 := 2 : 2 ( 2 :i, q), respectively; for q = 0.5 solid lines and for q G {0.01, 0.99} dotted lines. 


5 Small-sample efficiency - Misspecification of the random effects distribution 

An extensive simulation study is conducted (a) to gauge the small-sample efficiency of the ML and 
approximation in Kuss et al. [27]’s (hereafter KHS) methods, and (b) to investigate in detail the 
misspecification of the parametric margin or family of copulas of the random effects distribution. 

To simulate the data we have used the generation process in [42] to get heterogeneous study 
sizes; the simulation steps follow: 

1. Simulate the study size n from a shifted gamma distribution, i.e., n ~ sGamma(a = 1.2, /3 = 
0.01, lag = 30) and round off to the nearest integer. 

2. Simulate (ui, U 2 ) from a parametric family of copulas C(; r); r is converted to BVN, Frank 
and Clayton rotated by 90/270 dependence parameter 0 via the relations in (12), (13), and 
(14). 

3. Convert to beta realizations via Xj = F!-^(ttj,c)j, 7 j) or normal realizations via Xj = 
for j = 1, 2; for the latter convert to proportions via Xj = l~^{xj). 


11 






























4. Draw the number of diseased ni from a B{n, 0.43) distribution. 

5. Set 712 = n — Til, Uj = UjXj and then round Uj for j = 1, 2. 

We randomly generated B = 10^ samples of size N = 20, 50 from the Clayton rotated by 
270 degrees copula mixed model with beta margins. Table 2 contains the resultant biases, root 
mean square errors (RMSE), and standard deviations (SD), along with average theoretical variances 
scaled by N, for the MLEs under different copula choices and margins. The theoretical variances 
of the MEEs are obtained via the gradients and the Hessian computed numerically during the max¬ 
imization process. We also provide biases, RMSEs and SDs for the KHS estimates under the ‘true’ 
model, i.e., the Clayton rotated by 270 degrees copula mixed model with beta margins. 

Conclusions from the values in the table are the following: 

• ME with the the ‘true’ copula mixed model is highly efficient according to the simulated 
biases and standard deviations. 

• The MEEs of the meta-analytic parameters are slightly underestimated under copula mis- 
specification. That is, there is some downward bias for these parameters, especially if the 
“working” model is not close to Kullback-Eiebler distance with the “true” model, i.e., it is 
misspecified. Eor example in the table there is more bias for the Clayton rotated by 90 de¬ 
grees and Erank copulas since they have different tail dependence from the ‘true’ model, i.e., 
the rotated Clayton by 270 degrees. An interesting result is that the BVN copula performed 
rather well under misspecification. 

• The SDs are rather robust to the copula misspecification. 

• The meta-analytic MEEs and SDs are not robust to the margin misspecification, while the 
MEE of r and its SD is. 

• The KHS approximation method yields to biased univariate parameter estimates. 

• The efficiency of fhe KHS approximafion mefhod is low for fhe dependence paramefer r. The 
parameter r is subslanfially underesfimafed. 

The simulafion resulfs indicafe fhaf fhe KHS approximafion mefhod in [27] is an inefficient; 
hence flawed mefhod. This was expecfed, since fheorefically fhere are serious problems on mod¬ 
elling assumpfions under fhe case of heferogeneous sfudy sizes. If fhe number of frue positives 
and negafives do nof have a common supporf over differenl sfudies, fhen one cannof conclude fhaf 
fhere is a copula. To make our sfudy complete, we perform fheorefical calculafions, similarly fo 
[23, 35, 37], fo invesfigafe fhe accuracy of fhe approximafe copula likelihood mefhod in [27] for 
fhe special case of a consfanf size n of groups of diseased and healfhy people in fhe single studies 
and show whether or not this leads to consistent estimate of the parameters of the bivariate random 
effects distribution. As shown in the Appendix, the KHS method leads to asymptotic bias for both 
the univariate and copula parameters and hence there is no consistency. Also given the resultant 
substantial asymptotic downward bias for the dependence parameter, the approximation deterio¬ 
rates, and, hence cannot be used e.g., for prediction purposes via SROC curves. To this end, the 
KHS approximation method is not used in the sequel, since its inefficiency has been shown, and, if 
should be avoided for bivariate mefa-analysis of diagnostic fesf accuracy sfudies. 

The effecf of misspecifying fhe copula choice can be seen as minimal for bofh fhe univaiiafe 
paramefers and Kendall’s fau. However, nofe fhaf (a) fhe mefa-analyfic pai'amefers are a univariafe 
inference, and hence if is fhe univariate marginal disfribufion fhaf mailers and nof fhe lype of fhe 
copula, and, (b) as previously emphasized Kendall’s fau only accounfs for fhe dependence dom- 
inaled by fhe middle of fhe dafa (sensilivifies and specificities), and if is expecfed fo be similar 
amongsf differenl families of copulas. However, the tail dependence varies, as explained in Section 
4, and is a property to consider when choosing amongst different families of copulas, and, hence 
affects the shape of SROC curves, i.e., prediction. SROC will essentially show the effect of different 
model (random effect distribution) assumptions, since it is an inference that depends on the joint 
distribution. 
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Table 2: Small sample of sizes N = 20, 50 simulations ('10'* replications) from the Clayton rotated by 270 degrees copula mixed model with beta margins and resultant biases, root mean 
square errors ( RMSE) and standard deviations (SD), along with the square root of the average theoretical variances (\/V), scaled by N, for the MLEs under different copula choices and 
margins. We also provide biases, RMSEs and SDs for the KHS estimates under the ‘true’ model. 



Margin 

Copula 

TTl = 

- 20 

0.7 

N = 50 

II 

t= 

II 

0.9 

Af = 50 

71 = 
Al = 20 

0.2 

Af = 50 

72 = 
N = 2Q 

0.1 

N = 50 

II O 

t- 

II 

-0.5 

N = 50 

N Bias 

Beta 

BVN 

-0.02 

-0.09 

0.00 

-0.06 

-0.64 

-1.24 

-0.32 

-0.51 

-1.54 

-2.31 



Frank 

-0.23 

-0.72 

0.10 

0.27 

-0.59 

-1.08 

-0.39 

-0.72 

-2.01 

-3.77 



Clayton by 90 

-0.11 

-0.31 

-0.02 

-0.13 

-0.50 

-0.81 

-0.20 

-0.14 

-0.23 

2.74 



Clayton by 270 

-0.01 

-0.08 

0.03 

0.04 

-0.65 

-1.25 

-0.40 

-0.78 

-2.30 

-4.57 


Normal 

BVN 

0.71 

1.87 

0.63 

1.64 

- 

- 

- 

- 

-1.56 

-2.26 



Frank 

0.46 

1.11 

0.68 

1.79 

- 

- 

- 

- 

-1.95 

-3.52 



Clayton by 90 

0.62 

1.65 

0.64 

1.67 

- 

- 

- 

- 

-0.25 

2.66 



Clayton by 270 

0.70 

1.89 

0.62 

1.62 

- 

- 

- 

- 

-2.39 

-4.69 


KHS 

Clayton by 270 

1.99 

5.58 

-0.26 

-0.55 

-0.48 

-1.01 

-1.44 

-3.87 

7.42 

20.09 

A^SD 

Beta 

BVN 

0.94 

1.47 

0.44 

0.72 

1.01 

1.63 

0.71 

1.22 

3.31 

4.83 



Frank 

1.00 

1.58 

0.42 

0.68 

1.01 

1.62 

0.67 

1.13 

3.36 

4.84 



Clayton by 90 

0.97 

1.50 

0.46 

0.78 

1.10 

1.80 

0.87 

1.52 

4.83 

7.15 



Clayton by 270 

0.94 

1.47 

0.43 

0.70 

0.97 

1.56 

0.64 

1.07 

3.12 

4.54 


Normal 

BVN 

1.06 

1.66 

0.38 

0.59 

4.77 

7.35 

5.02 

7.93 

3.35 

4.76 



Frank 

1.11 

1.76 

0.38 

0.59 

4.70 

7.24 

5.04 

7.99 

3.40 

4.92 



Clayton by 90 

1.11 

1.72 

0.38 

0.60 

5.19 

8.13 

5.52 

8.88 

4.76 

6.74 



Clayton by 270 

1.06 

1.66 

0.38 

0.59 

4.52 

6.86 

4.97 

7.77 

3.18 

4.62 


KHS 

Clayton by 270 

1.17 

1.95 

0.56 

0.77 

1.04 

1.68 

0.76 

0.72 

1.82 

1.48 

nVv 

Beta 

BVN 

0.76 

1.25 

0.40 

0.66 

0.79 

1.29 

0.61 

1.02 

2.67 

3.86 



Frank 

0.77 

1.27 

0.37 

0.59 

0.82 

1.33 

0.58 

0.95 

2.62 

3.78 



Clayton by 90 

0.77 

1.25 

0.39 

0.65 

0.83 

1.34 

0.64 

1.10 

3.00 

4.11 



Clayton by 270 

0.75 

1.22 

0.38 

0.60 

0.77 

1.25 

0.56 

0.90 

2.56 

3.72 


Normal 

BVN 

0.87 

1.45 

0.32 

0.53 

3.70 

5.82 

4.58 

7.35 

2.90 

3.84 



Frank 

0.87 

1.45 

0.31 

0.50 

3.74 

5.92 

4.58 

7.31 

2.56 

3.69 



Clayton by 90 

0.89 

1.48 

0.33 

0.54 

3.87 

6.15 

4.57 

7.44 

2.90 

3.94 



Clayton by 270 

0.84 

1.36 

0.31 

0.49 

3.52 

5.47 

4.46 

6.96 

2.45 

3.48 

AT RMSE 

Beta 

BVN 

0.94 

1.47 

0.44 

0.72 

1.19 

2.05 

0.78 

1.32 

3.65 

5.35 



Frank 

1.02 

1.74 

0.43 

0.73 

1.17 

1.95 

0.77 

1.34 

3.92 

6.13 



Clayton by 90 

0.97 

1.53 

0.46 

0.80 

1.21 

1.97 

0.89 

1.53 

4.84 

7.65 



Clayton by 270 

0.94 

1.47 

0.43 

0.70 

1.17 

2.00 

0.75 

1.33 

3.87 

6.44 


Normal 

BVN 

1.27 

2.50 

0.74 

1.74 

- 

- 

- 

- 

3.69 

5.27 



Frank 

1.20 

2.08 

0.78 

1.88 

- 

- 

- 

- 

3.92 

6.05 



Clayton by 90 

1.27 

2.39 

0.75 

1.77 

- 

- 

- 

- 

4.77 

7.24 



Clayton by 270 

1.27 

2.52 

0.73 

1.73 

- 

- 

- 

- 

3.98 

6.58 


KHS 

Clayton by 270 

2.31 

5.91 

0.62 

0.94 

1.15 

1.96 

1.63 

3.93 

7.64 

20.14 








6 Vuong’s test for model comparison 

In this section we provide a methodology for the comparison of non-nested models. It would be 
used as a tool to show if the copula mixed model provides better fit than the standard GLMM. We 
will call a test proposed by Vuong [53]. The Vuong’s test is the sample version of the difference in 
Kullback-Leibler divergence between two models and can be used to differentiate two parametric 
models which could be non-nested.This test has been used extensively in the copula literature to 
compare copula models, see e.g., [2, 3, 25] 

Assume that we have Models 1 and 2 with parametric densities and respectively. We 
can compare 


— N 1 [log /^(Vi, Va)] - [log V 2 ; 

i 


and 

= N-^[Y^{Ef^[logf^{Yi,Y2)] - Ej^[\ogf^^\YuY2-,e^^^)]} 

i 


where , 0^^^ are the parameters in Models 1 and 2 respectively that lead to the closest Kullback- 
Leibler divergence to the true equivalently they are the limits in probability of the MLEs based 
on models 1 and 2 respectively. Model 1 is closer to the true /^, i.e., is the better fitting model if 
A = — A 2 y'i' < 0, and Model 2 is the better fitting model if A > 0. The sample version of 

A with MLEs 0^^\ 0^^^ is 

N 

D = J2Di/N, 

i=l 


where Di 


log 


/(i)(yi,y2;0‘'^) 


Vuong [53] has shown that asymptotically 


VnD/s ~ A(0,1), 

where ^ E*=i(A - D?- 


1 Illustrations 

We illustrate the use of the copula mixed model by re-analysing the data of two published meta¬ 
analysis [48, 13]. These data have been frequently used as an example for methodological papers 
on meta-analysis of diagnostic accuracy studies [43, 4, 46, 44, 16, 27]. 

We fit the copula mixed model for all different choices of parametric families of copulas and 
margins. To make it easier to compare strengths of dependence, we convert the copula parameters to 
Kendall’s r’s via the relations in (12), (13), and (14) for BVN, Erank and rotated Clayton copulas, 
respectively. Since the number of parameters is the same between the models, we use the log- 
likelihood at estimates as a rough diagnostic measure for goodness of fit between the models. We 
further compute the Vuong’s tests with Model 1 being the BVN copula mixed model with normal 
margins, i.e., the standard GLMM, to reveal if any other copula mixed model provides better fit than 
the standard GLMM. 

Einally, we demonstrate SROC curves and summary operating points (a pair of average sensi¬ 
tivity and specificity) with a confidence region and a predicfive region as deduced in Secfion 4.1. 


7.1 The telomerase and computed tomography data 

In Glas et al. [13] the telomerase marker for the diagnosis of bladder cancer is evaluated using 10 
studies. The size in each study ranges from 35 to 195. The interest was to define if this non-invasive 
and cheap marker could replace the standard of cystoscopy or histopathology. Riley et al. [44] 
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applied the GLMM with different starting values and all produced a between-study correlation esti¬ 
mate of — 1 but with different meta-analytic parameter point estimates and standard errors. Clearly 
at this example, it is not possible to estimate the correlation between the logit sensitivity and speci¬ 
ficity, and the maximum likelihood estimator should truncate the correlation to the left boundary of 
its parameter space, i.e., —1. In [27] it is acknowledged that the copula model for observed discrete 
variables which have beta-binomial margins yields sensible estimates for the dependence parameter 
and its standard error. This result is in error and due to the fact that the KHS method underestimates 
the dependence parameter as emphasized in Section 5 and shown in the Appendix for p = —1. 

Fitting the copula mixed model for all different choices of parametric families of copulas and 
margins, the resultant estimate of the dependence parameter was close to the boundary of its param¬ 
eter space. If the dependence parameter estimate is so large (on absolute value), the copula should 
be set to countermonotonic (Frechet lower bound), and, then optimize over the remaining (univari¬ 
ate) parameters. With other words there exists negative perfect dependence, and thus there is only 
one copula: the countermonotonic copula. This is a limiting case for all the parametric families 
of copulas, listed in Section 4, when the dependence parameter is fixed fo fhe leff boundary of ifs 
paramefer space. 

This was also fhe case for fhe analysis of fhe dafa on 17 sfudies of compufed fomography (CT) 
for fhe diagnosis of lymph node mefasfasis in women wifh cervical cancer, one of fhree imaging 
fechniques in fhe mefa-analysis in [48]. The size in each sfudy ranges from 20 fo 253. Diagnosis of 
mefasfafic disease by CT relies on nodal enlargemenf. 


Table 3: Maximised log-likelihoods, estimates and standard errors (SE), along with the Vuong’s statistics and p-values 
for the telomerase and computed tomography data. 


Telomerase 


Computed Tomography 


Normal 

margins 

Beta 

margins 

Normal margins 

Beta 

margins 


Est. 

SE 


Est. 

SE 


Est. 

SE 


Est. 

SE 

TTl 

0.77 

0.03 

TTl 

0.76 

0.03 

TTl 

0.46 

0.07 

TTl 

0.46 

0.06 

7r2 

0.91 

0.05 

7r2 

0.81 

0.06 

7V2 

0.93 

0.01 

7V2 

0.92 

0.01 

CJ-1 

0.43 

0.13 

71 

0.03 

0.02 

(Jl 

1.00 

0.27 

71 

0.17 

0.07 

0-2 

1.83 

0.40 

72 

0.28 

0.10 

(72 

0.60 

0.23 

72 

0.02 

0.02 

logL 

-50.37 

logL 

-51.14 

logL 

-69.37 

logL 

-69.58 



Vuong’s test 





Vuonj 

'’s test 



s/ND/s 



VND/s 

-1.580 





-1.416 

p-value 



p-value 

0.114 





0.157 


Table 3 gives fhe esfimafed univariafe paramefers, sfandard errors, and log-likelihoods for bofh 
normal and befa margins for bofh dafasefs. For felomerase dafa, bofh models agree on fhe esfi- 
mafed sensifivify ffi buf fhe esfimafe of specificify tt 2 is larger under fhe sfandard CLMM. The 
log-likelihood is —50.37 for normal margins and —51.14 for befa margins, and fhus a normal mar¬ 
gin seems fo be a heller hi for fhe dafa. Furthermore, according fo fhe Vuong’s lesf fhe copula 
mixed model wifh normal margins (i.e., fhe sfandard GLMM) provides marginally heller fif {p- 
value= 0.114). For computed tomography dafa, bofh models agree on fhe esfimafed sensifivify tti 
and specificify 7 ( 2 . The log-likelihood is —69.37 for normal margins and —69.58 for befa margins, 
and fhus a normal margin seems to be a belter fif for fhe dafa. However, according to fhe Vuong’s lesf 
fhe copula mixed model wifh normal margins (i.e., fhe sfandard GLMM) does nof provide slalislical 
significanl heifer fif (p-value= 0.157). 

Finally, figure 3 also shows fhe SROC curves for bofh dafasefs and fhe visual fif is consisfenf 
wifh fhe model filling and comparison in Table 3. Note in passing since we are dealing wifh fhe 
counfermonofonic copula all fhe quantile regression curves almosf coincide, and hence we jusf 
depicl one median regression curve for each model. 
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Figure 3: SROC curves from the countermonotonic copula representation of the random effects distribution with normal 
margins (black line) and beta (red line) margins for the telomerase and computed tomography data. 


7.2 The lymphangiography data 

In this section we apply the copula mixed models to data on 17 studies of lymphangiography for the 
diagnosis of lymph node metastasis in women with cervical cancer, one of three imaging techniques 
in the meta-analysis in [48]. The size in each study ranges from 21 to 300. Diagnosis of metastatic 
disease by lymphangiography is based on the presence of nodal-filling defects. 


Table 4: Maximised log-likelihoods, estimates and standard errors (SE), along with the Vuong’s statistics and p-values 
for the lymphangiography data. 


Normal margins 



BVN 

Estimate 

SE 

Frank 

Estimate 

SE 

Clayton by 180 
Estimate SE 

Clayton by 270 
Estimate SE 

TTl 

0.67 

0.03 

0.68 

0.03 

0.67 

0.03 

0.67 

0.03 

7r2 

0.84 

0.03 

0.84 

0.03 

0.84 

0.03 

0.84 

0.04 

CTl 

0.35 

0.19 

0.36 

0.18 

0.34 

0.18 

0.34 

0.19 

0'2 

0.91 

0.22 

0.91 

0.22 

0.91 

0.22 

0.90 

0.22 

T 

0.16 

0.29 

0.18 

0.28 

0.14 

0.21 

0.19 

0.29 

log 7/ 

-91.38 


-91.32 


-91.32 


-91.15 


Vuong’s test 

s/ND/s 

- 


0.523 


0.274 


1.280 


p-value 

- 


0.601 


0.784 


0.201 


Beta margins 


BVN 


Frank 


Clayton by 180 

Clayton by 270 


Estimate 

SE 

Estimate 

SE 

Estimate 

SE 

Estimate 

SE 

TTl 

0.67 

0.03 

0.67 

0.03 

0.67 

0.03 

0.67 

0.03 

7r2 

0.81 

0.03 

0.81 

0.03 

0.81 

0.03 

0.81 

0.03 

7i 

0.03 

0.03 

0.03 

0.03 

0.02 

0.03 

0.02 

0.03 

72 

0.09 

0.04 

0.09 

0.04 

0.10 

0.04 

0.09 

0.04 

T 

0.15 

0.30 

0.18 

0.32 

0.16 

0.40 

0.19 

0.28 

log 7/ 

-90.67 


-90.61 


-90.60 


-90.44 


Vuong’s test 

yjND/s 

1.668 


1.798 


1.877 


2.248 


p-value 

0.095 


0.072 


0.061 


0.025 



In Table 4 we report the resulting maximized log-likelihoods, estimates, and standard errors of 
the copula mixed models with different choices of parametric families of copulas and margins. All 
models agree on the estimated sensitivity tti, but the estimate ft 2 of specificity is smaller when beta 
margins are assumed. The log-likelihoods show that a copula mixed model with rotated by 270 
degrees Clayton copula and beta margins provides the best fit. It is revealed that a copula mixed 
model with the sensitivity and specificity on the original scale provides better fit than the GLMM, 
which models the sensitivity and specificity on a transformed scale. The improvement over the 
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Figure 4: Contour plots and quantile regression curves from the copula representation of the random effects distribution 
with normal margins and BVN, Frank, and Clayton by 180 and 270 copulas for the lymphangiography data. Red and 
green lines represent the quantile regression curves xi := '£ 1 ( 0 : 2 , q) and X2 ■= X2{xi, q), respectively; for q — 0.5 solid 
lines and for q G {0.01, 0.99} dotted lines. 


GLMM is small in terms of the likelihood principle, but for the Vuong’s statistic there is enough 
improvement to get a statistical significant difference (p-value= 0.025). 

Figures 4 and 5 show the fitted SROC curves along with their confidence and prediction regions 
for the copula mixed models with normal and beta margins, respectively. Note that the predictive 
regions cover a greater range of specificity rather than sensitivity. 


7.3 The magnetic resonance imaging data 

In this section we apply the copula mixed models to data on 10 studies of magnetic resonance 
imaging for the diagnosis of lymph node metastasis in women with cervical cancer, the last imaging 
technique in the meta-analysis in [48]. The size in each study ranges from 20 to 272. Diagnosis of 
metastatic disease by lymphangiography relies on nodal enlargement. 

In Table 5 we report the resulting maximized log-likelihoods, estimates, and standard errors of 
the copula mixed models with different choices of parametric families of copulas and margins. All 
models roughly agree on the estimated sensitivity -tti and specificity 7 ( 2 , but both are slightly smaller 
when beta margins are assumed. The log-likelihoods show that a rotated by 270 degrees Clayton 
copula mixed model with normal or beta margins provides the best fit. Although, the rotated by 270 
degrees Clayton copula mixed model provides better fit than the GLMM, the difference, according 
to Vuong’s test, is not statistical significant (p-value=0.156). 

Figures 6 and 7 show the fitted SROC curves along with their confidence and predicfion regions 
for fhe copula mixed models wifh normal and befa margins, respecfively. Nofe fhaf fhe predicfive 
regions cover a greafer range of sensifivify rafher fhan specificily. 
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Figure 5: Contour plots and quantile regression curves from the copula representation of the random effects distribution 
with beta margins and BVN, Frank, and Clayton by 180 and 270 copulas for the lymphangiography data. Red and green 
lines represent the quantile regression curves xi := xi{x 2 , q) and X 2 ■= X 2 {xi, q), respectively; for q = 0.5 solid lines 
and for q £ {0.01, 0.99} dotted lines. 

Table 5: Maximised log-likelihoods, estimates and standard errors (SE), along with the Vuong's statistics and p-values 
for the magnetic resonance imaging data. 


Normal margins 



BVN 

Estimate 

SE 

Frank 

Estimate 

SE 

Clayton by 90 
Estimate SE 

Clayton by 270 
Estimate SE 

TTl 

0.55 

0.11 

0.54 

0.10 

0.54 

0.11 

0.55 

0.10 

7r2 

0.95 

0.02 

0.96 

0.02 

0.95 

0.02 

0.96 

0.02 

(Tl 

1.16 

0.39 

1.14 

0.38 

1.21 

0.41 

1.13 

0.37 

0-2 

0.87 

0.34 

0.83 

0.32 

0.85 

0.34 

0.87 

0.32 

T 

-0.51 

0.29 

-0.47 

0.28 

-0.48 

0.33 

-0.49 

0.26 

logL 

-46.26 


-46.35 


-46.72 


-45.90 


Vuong’s test 

yj ND/S 

- 


-0.815 


-2.175 


1.419 


p-value 

- 


0.415 


0.030 


0.156 


Beta margins 


BVN 


Frank 


Clayton by 90 

Clayton by 270 


Estimate 

SE 

Estimate 

SE 

Estimate 

SE 

Estimate 

SE 

TTl 

0.54 

0.08 

0.53 

0.08 

0.53 

0.08 

0.54 

0.08 

7r2 

0.94 

0.02 

0.94 

0.02 

0.94 

0.02 

0.94 

0.02 

71 

0.21 

0.10 

0.21 

0.10 

0.22 

0.10 

0.21 

0.09 

72 

0.04 

0.03 

0.03 

0.02 

0.03 

0.03 

0.04 

0.02 

T 

-0.53 

0.28 

-0.47 

0.28 

-0.50 

0.33 

-0.50 

0.25 

logL 

-46.27 


-46.39 


-46.75 


-45.86 


Vuong’s test 

s/ND/s 

-0.014 


-0.422 


-1.326 


0.935 


p-value 

0.989 


0.673 


0.185 


0.350 
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Figure 6: Contour plots and quantile regression curves from the copula representation of the random effects distribution 
with normal margins and BVN, Frank, and Clayton by 90 and 270 copulas for the magnetic resonance imaging data. Red 
and green lines represent the quantile regression curves xi := iEi(a: 2 , q) and X 2 ■= X 2 {xi, q), respectively: for q = 0.5 
solid lines and for q G {0.01, 0.99} dotted lines. 


8 Discussion 

We have proposed a copula mixed model for bivariate meta-analysis of diagnostic test accuracy 
studies. This is the most general meta-analytic model, with univariate parameters separated from 
dependence parameters. Our general model includes the GLMM as a special case and can provide 
an improvement over the latter based on log-likelihood and Vuong’s [53] statistic, and thus can 
provide a better statistical inference for the SROC. This improvement relies on the fact that the 
random effects distribution is expressed via copulas which allow for flexible dependence modelling, 
different from assuming simple linear correlation structures, normality and tail independence, which 
makes them well suited to the aforementioned application area. 

Building on the basic model proposed in this paper, there are several extensions that can be 
implemented. The copula mixed model can also easily be extended in any context where clinical 
trials or observational studies report more than a single outcome and to inclusion of covariates. 
However, larger sample sizes will be required to estimate the effect of covariates in bivariate meta¬ 
regression, where the underlying treatment effects depend on covariates. This is typical in the 
univariate meta-regression [20]. 

Another direction of future research is to extend our copula-based meta-analytic model to the 
d-variate {d > 2) case. There are many simple bivariate copula families, but generally their mul¬ 
tivariate extensions have limited dependence structures. However, in recent years, a popular and 
useful approach is the vine pair-copula construction, see e.g., [26, 25], which is based on d{d— 1)/2 
bivariate copulas. Some studies also may not report all d outcomes. In such cases our model can be 
extended for missing data via pattern mixture models. Pattern mixture models are studied in [50] 
for copulas and in [31] for pairwise and network meta-analysis. 
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Figure 7: Contour plots and quantile regression curves from the copula representation of the random effects distribution 
with beta margins and BVN, Frank, and Clayton by 90 and 270 copulas for the magnetic resonance imaging data. Red 
and green lines represent the quantile regression curves xi := iEi(a: 2 , q) and X 2 ■= X 2 {xi, q), respectively; for q — 0.5 
solid lines and for q G {0.01, 0.99} dotted lines. 


Software 

A contributed R package CopulaREMADA [36] has functions to implement the copula mixed model 
for meta-analysis of diagnostic test accuracy studies and produce SROC curves and summary oper¬ 
ating points (a pair of average sensitivity and specificity) with a confidence region and a predictive 
region. All the analyses presented in Section 7 are given as code examples in the package. The 
R package VGAM [54] and specifically the functions pbetabinom and dbetabinom have been 
used to implement the marginal distributions for the KHS approximation method in [27]. 

Appendix 

We study the asymptotics of the KHS approximation method in [27], and we assess the accuracy 
based on the limit (as the number of clusters increases to infinity) of the maximum KHS likelihood 
estimate (KHSMLE). By varying factors such as the marginal and copula parameters we demon¬ 
strate patterns in the asymptotic bias of the KHSMLE, and assess the performance of KHS . We 
will compute these limiting KHSMLE in a variety of situations to show clearly if the KHS method 
is good. By using this limit, we show whether or not this leads to consistent estimate of the param¬ 
eters of the bivariate random effects distribution; hence prove whether the KHS approach is valid 
or not. Eor the cases where we compute the probability limit, we will take a constant size n of 
groups of diseased and healthy people in the single studies that increases. Eor ease of exposition, 
we also consider the case that the univariate marginal parameters are common to different univariate 
margins. 

Let the T distinct cases for the discrete response be denoted as 
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In a random sample of size N, let the corresponding frequencies be denoted as , ■ ■ ■, . Let 

be the limit in probability of /N as —)• oo. 


Table Al: Limiting KHSMLEfor a BVN copula mixed model with beta margins. 


p 

P 

ITFrS - 

TT 

TT 

ITH^ 

7 

7 

ITFT^ 


n = 20 

n = 100 


n = 20 

n = 100 


n = 20 

n = 100 

-0.2 

-0.047 

-0.163 

0.7 

0.701 

0.701 

0.05 

0.050 

0.050 

-0.5 

-0.160 

-0.412 

0.7 

0.705 

0.702 

0.05 

0.049 

0.050 

-0.8 

-0.277 

-0.664 

0.7 

0.709 

0.703 

0.05 

0.046 

0.050 

-1 

-0.356 

-0.825 

0.7 

0.711 

0.704 

0.05 

0.044 

0.049 

-0.2 

-0.051 

-0.174 

0.7 

0.702 

0.701 

0.1 

0.099 

0.100 

-0.5 

-0.164 

-0.438 

0.7 

0.708 

0.703 

0.1 

0.095 

0.099 

-0.8 

-0.289 

-0.708 

0.7 

0.714 

0.705 

0.1 

0.085 

0.096 

-1 

-0.381 

-0.885 

0.7 

0.718 

0.706 

0.1 

0.075 

0.089 

-0.2 

-0.030 

-0.122 

0.7 

0.703 

0.702 

0.2 

0.199 

0.198 

-0.5 

-0.129 

-0.320 

0.7 

0.715 

0.707 

0.2 

0.187 

0.187 

-0.8 

-0.242 

-0.558 

0.7 

0.728 

0.714 

0.2 

0.162 

0.161 

-1 

-0.338 

-0.788 

0.7 

0.738 

0.722 

0.2 

0.133 

0.116 

-0.2 

-0.007 

-0.155 

0.8 

0.800 

0.801 

0.05 

0.050 

0.050 

-0.5 

-0.076 

-0.396 

0.8 

0.804 

0.802 

0.05 

0.049 

0.049 

-0.8 

-0.148 

-0.642 

0.8 

0.808 

0.804 

0.05 

0.045 

0.048 

-1 

-0.191 

-0.802 

0.8 

0.811 

0.805 

0.05 

0.041 

0.045 

-0.2 

-0.005 

-0.130 

0.8 

0.800 

0.801 

0.1 

0.100 

0.099 

-0.5 

-0.084 

-0.339 

0.8 

0.808 

0.804 

0.1 

0.095 

0.094 

-0.8 

-0.168 

-0.575 

0.8 

0.817 

0.808 

0.1 

0.084 

0.082 

-1 

-0.230 

-0.765 

0.8 

0.823 

0.813 

0.1 

0.073 

0.064 

-0.2 

0.026 

-0.062 

0.8 

0.795 

0.803 

0.2 

0.202 

0.196 

-0.5 

-0.064 

-0.187 

0.8 

0.814 

0.812 

0.2 

0.188 

0.182 

-0.8 

-0.165 

-0.348 

0.8 

0.837 

0.825 

0.2 

0.156 

0.148 

-1 

-0.248 

-0.535 

0.8 

0.853 

0.837 

0.2 

0.122 

0.104 

-0.2 

0.084 

-0.075 

0.9 

0.890 

0.901 

0.05 

0.052 

0.049 

-0.5 

0.031 

-0.214 

0.9 

0.896 

0.904 

0.05 

0.051 

0.046 

-0.8 

-0.025 

-0.377 

0.9 

0.903 

0.907 

0.05 

0.048 

0.039 

-1 

-0.062 

-0.515 

0.9 

0.908 

0.910 

0.05 

0.045 

0.031 

-0.2 

0.114 

-0.035 

0.9 

0.878 

0.902 

0.1 

0.110 

0.098 

-0.5 

0.045 

-0.141 

0.9 

0.891 

0.908 

0.1 

0.105 

0.089 

-0.8 

-0.034 

-0.269 

0.9 

0.908 

0.916 

0.1 

0.092 

0.071 

-1 

-0.092 

-0.396 

0.9 

0.920 

0.923 

0.1 

0.078 

0.051 

-0.2 

0.182 

0.032 

0.9 

0.837 

0.893 

0.2 

0.241 

0.193 

-0.5 

0.127 

-0.054 

0.9 

0.858 

0.911 

0.2 

0.225 

0.168 

-0.8 

0.050 

-0.179 

0.9 

0.894 

0.936 

0.2 

0.178 

0.115 

-1 

-0.032 

-0.324 

0.9 

0.926 

0.954 

0.2 

0.126 

0.061 


For the KHS log-likelihood we have the limit, 


N ^£(^,7,0)log 


t=l 


Z 

'[H{yf^-,n,TT,'y),H{y^2^-,n,Tr,j);e'j n,vr,7) 


i=i 


( 15 ) 


The limit of the KHSMLE (as N —>• 00) is the maximum of ( 15 ); we denote this limit as ^ 

The in ( 15 ) are the model based probabilities and are computed to at least five significant digits 
using Gauss-Legendre quadrature [ 52 ] with a sufficient number of quadrature points as described in 
Subsection 3 . 4 . For the log-likelihood in ( 9 ), we have the limit, 



c(ni, U2; 9 )duidu 2 - 


( 16 ) 


The limit of the MLE (as N —)• 00) is the maximum of ( 16 ); we denote this limit as (vr, 7, f). 
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Representative results are shown in Table A1 for a BVN copula mixed model with beta margins, 
with MLE results omitted because they were identical with the true values up to four or five decimal 
places. Therefore, our method leads to unbiased estimating equations. Regarding the KHS method, 
conclusions from the values in the table and other computations that we have done are that for the 
KHS method there is asymptotic bias (decreases as n increases) for the univariate parameters vr and 
7 as TT, 7 and p increase, and substantial asymptotic downward bias for the dependence parameter 
/?; note that this slightly decreases as n increases. 
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